library(xtable)

quartz(width=9,height=6)

par(mfrow = c(2, 3), pty = "s",family="Gill Sans MT",font.main=1,cex=.9)#,cex.lab=0.9,cex=1.1,cex.axis=0.9)

#### Right

load("~/Dropbox/complieropt/Polmeth Archive/Right.rdata")

simoutput[8,] == 500 -> isN
t500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 1000 -> isN
t1000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 2500 -> isN
t2500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 5000 -> isN
t5000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 10000 -> isN
t10000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 25000 -> isN
t25000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))

plot((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,4],x=c(500,1000,2500,5000,10000,25000),type="p",ylim=c(0,3),pch=1,xlab="N",ylab="RMSE",main="Right Model")

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,3],x=c(500,1000,2500,5000,10000,25000),pch=2)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,2],x=c(500,1000,2500,5000,10000,25000),pch=3)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,1],x=c(500,1000,2500,5000,10000,25000),pch=4)

xtable(rbind(t500,t1000,t2500,t5000,t10000,t25000))

load("~/Dropbox/complieropt/Polmeth Archive/Final Code/Noise.rdata")

simoutput[8,] == 500 -> isN
t500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 1000 -> isN
t1000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 2500 -> isN
t2500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 5000 -> isN
t5000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 10000 -> isN
t10000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 25000 -> isN
t25000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))

xtable(rbind(t500,t1000,t2500,t5000,t10000,t25000))

plot((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,4],x=c(500,1000,2500,5000,10000,25000),type="p",ylim=c(0,3),pch=1,xlab="N",ylab="RMSE",main="Noise (Heterogeneous)")

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,3],x=c(500,1000,2500,5000,10000,25000),pch=2)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,2],x=c(500,1000,2500,5000,10000,25000),pch=3)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,1],x=c(500,1000,2500,5000,10000,25000),pch=4)

simoutput[8,] == 500 -> isN
t500 <- c(colMeans(apply(simoutput[9:12,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 1000 -> isN
t1000 <- c(colMeans(apply(simoutput[9:12,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 2500 -> isN
t2500 <- c(colMeans(apply(simoutput[9:12,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 5000 -> isN
t5000 <- c(colMeans(apply(simoutput[9:12,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 10000 -> isN
t10000 <- c(colMeans(apply(simoutput[9:12,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 25000 -> isN
t25000 <- c(colMeans(apply(simoutput[9:12,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))

xtable(rbind(t500,t1000,t2500,t5000,t10000,t25000))

plot((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,4],x=c(500,1000,2500,5000,10000,25000),type="p",ylim=c(0,3),pch=1,xlab="N",ylab="RMSE",main="Noise (Homogeneous)")

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,3],x=c(500,1000,2500,5000,10000,25000),pch=2)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,2],x=c(500,1000,2500,5000,10000,25000),pch=3)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,1],x=c(500,1000,2500,5000,10000,25000),pch=4)

load("~/Dropbox/complieropt/Polmeth Archive/Final Code/Wrong Model.rdata")

simoutput[8,] == 500 -> isN
t500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 1000 -> isN
t1000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 2500 -> isN
t2500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 5000 -> isN
t5000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 10000 -> isN
t10000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 25000 -> isN
t25000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))

xtable(rbind(t500,t1000,t2500,t5000,t10000,t25000))

plot((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,4],x=c(500,1000,2500,5000,10000,25000),type="p",ylim=c(0,3),pch=1,xlab="N",ylab="RMSE",main="Wrong Model")

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,3],x=c(500,1000,2500,5000,10000,25000),pch=2)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,2],x=c(500,1000,2500,5000,10000,25000),pch=3)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,1],x=c(500,1000,2500,5000,10000,25000),pch=4)

load("~/Dropbox/complieropt/Polmeth Archive/Final Code/Only X1.rdata")

simoutput[8,] == 500 -> isN
t500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 1000 -> isN
t1000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 2500 -> isN
t2500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 5000 -> isN
t5000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 10000 -> isN
t10000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 25000 -> isN
t25000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))

xtable(rbind(t500,t1000,t2500,t5000,t10000,t25000))


plot((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,4],x=c(500,1000,2500,5000,10000,25000),type="p",ylim=c(0,3),pch=1,xlab="N",ylab="RMSE",main="Only X1")

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,3],x=c(500,1000,2500,5000,10000,25000),pch=2)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,2],x=c(500,1000,2500,5000,10000,25000),pch=3)

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,1],x=c(500,1000,2500,5000,10000,25000),pch=4)


plot.new()

legend(x=0,y=1,legend=c("OLS (Z)","OLS (D)", "2SLS", "ICSW"),pch=c(4,3,2,1),cex=1.2) # double check


#### Permutation Tests

quartz(width=6,height=6)

par(mfrow = c(1, 1), pty = "s",family="Gill Sans MT",font.main=1)#,cex.lab=0.9,cex=1.1,cex.axis=0.9)

load("~/Dropbox/complieropt/Polmeth Archive/Final Code/Right.rdata")

simoutput[8,] == 500 -> isN
t500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 1000 -> isN
t1000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 2500 -> isN
t2500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 5000 -> isN
t5000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 10000 -> isN
t10000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 25000 -> isN
t25000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))

plot((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,5],x=c(500,1000,2500,5000,10000,25000),type="p",ylim=c(0,0.75),pch=1,xlab="N",ylab="Average p",main="Permutation Test Performance")

load("~/Dropbox/complieropt/Polmeth Archive/Final Code/Noise.rdata")

simoutput[8,] == 500 -> isN
t500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 1000 -> isN
t1000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 2500 -> isN
t2500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 5000 -> isN
t5000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 10000 -> isN
t10000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 25000 -> isN
t25000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))

points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,5],x=c(500,1000,2500,5000,10000,25000),type="p",ylim=c(0,1),pch=2,xlab="N",ylab="Average p",main="Permutation Test Performance")

load("~/Dropbox/complieropt/Polmeth Archive/Final Code/Wrong Model.rdata")

simoutput[8,] == 500 -> isN
t500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 1000 -> isN
t1000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 2500 -> isN
t2500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 5000 -> isN
t5000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 10000 -> isN
t10000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 25000 -> isN
t25000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))


points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,5],x=c(500,1000,2500,5000,10000,25000),type="p",ylim=c(0,1),pch=3,xlab="N",ylab="Average p",main="Permutation Test Performance")

load("~/Dropbox/complieropt/Polmeth Archive/Final Code/Only X1.rdata")

simoutput[8,] == 500 -> isN
t500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 1000 -> isN
t1000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 2500 -> isN
t2500 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 5000 -> isN
t5000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 10000 -> isN
t10000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))
simoutput[8,] == 25000 -> isN
t25000 <- c(colMeans(apply(simoutput[2:5,isN],1,function(x) (x - simoutput[1,isN])^2))^.5,mean(simoutput[7,isN]))


points((rbind(t500,t1000,t2500,t5000,t10000,t25000))[,5],x=c(500,1000,2500,5000,10000,25000),type="p",ylim=c(0,1),pch=4,xlab="N",ylab="Average p",main="Permutation Test Performance")


legend(x=16500,y=.75,legend=c("Correct","Noise", "Wrong Model", "Only X1"),pch=c(1,2,3,4),cex=1) # double check

